#!/usr/bin/env julia
"""
Test Script for 3D Lid-Driven Cavity Solver

This script tests the 3D Navier-Stokes solver implementation with a small problem
to verify that all the mathematical formulation, algorithms, and integrations work correctly.

Run from the project root directory:
```bash
cd LidDrivenCavity.jl
julia examples/test_3d_solver.jl
```
"""

using Pkg
Pkg.activate(".")

# Load the package
using LidDrivenCavity
using Printf

function test_3d_solver_basic()
    println("="^60)
    println("Testing 3D Lid-Driven Cavity Solver")
    println("="^60)
    
    # Small problem for testing
    options = LidDrivenOptions(
        n = 8,           # Very small grid (10x10x10 total nodes)
        Re = 100.0,      # Moderate Reynolds number
        max_steps = 50,  # Limited steps for testing
        tol = 1e-4,      # Relaxed tolerance for quick test
        solver = :julia, # Use built-in solver (most compatible)
        verbose = true,  # Show detailed output
        U_lid = 1.0      # Unit lid velocity
    )
    
    println("\n📋 Test Configuration:")
    println("   Grid size: $(options.n+2)³ = $((options.n+2)^3) total nodes")  
    println("   L-shape removes: ~25% of nodes")
    println("   Fluid nodes: ~$((options.n+2)^3 * 0.75) nodes")
    println("   Reynolds number: $(options.Re)")
    println("   Max iterations: $(options.max_steps)")
    println("   Tolerance: $(options.tol)")
    
    # Run the solver
    println("\n🚀 Starting 3D solver test...")
    
    try
        result = solve_lid_driven_3d(options)
        
        # Check basic result properties
        println("\n✅ Solver completed successfully!")
        println("   Converged: $(result.converged)")
        println("   Iterations: $(result.iterations)")  
        println("   Final residual: $(@sprintf("%.3e", result.residual_norm))")
        println("   Solve time: $(@sprintf("%.2f", result.solve_time))s")
        
        # Analyze solution properties
        analyze_3d_solution(result)
        
        return true, result
        
    catch e
        println("\n❌ Test failed with error:")
        println("   Error: $e")
        
        if isa(e, BoundsError)
            println("   This might be an indexing error in the 3D implementation.")
        elseif isa(e, DimensionMismatch) 
            println("   This might be a dimension mismatch in array operations.")
        elseif isa(e, MethodError)
            println("   This might be a missing function or incorrect arguments.")
        end
        
        println("\n🔍 Stack trace:")
        Base.show_backtrace(stdout, catch_backtrace())
        
        return false, nothing
    end
end

function analyze_3d_solution(result::LidDrivenResult)
    println("\n📊 Solution Analysis:")
    
    # Basic field statistics
    u_stats = (min=minimum(result.u), max=maximum(result.u), mean=mean(result.u))
    v_stats = (min=minimum(result.v), max=maximum(result.v), mean=mean(result.v))  
    w_stats = (min=minimum(result.w), max=maximum(result.w), mean=mean(result.w))
    
    println("   u-velocity: min=$(@sprintf("%.3f", u_stats.min)), max=$(@sprintf("%.3f", u_stats.max)), mean=$(@sprintf("%.3f", u_stats.mean))")
    println("   v-velocity: min=$(@sprintf("%.3f", v_stats.min)), max=$(@sprintf("%.3f", v_stats.max)), mean=$(@sprintf("%.3f", v_stats.mean))")
    println("   w-velocity: min=$(@sprintf("%.3f", w_stats.min)), max=$(@sprintf("%.3f", w_stats.max)), mean=$(@sprintf("%.3f", w_stats.mean))")
    
    # Check for expected lid-driven cavity behavior
    println("\n🔬 Physical validation:")
    
    # Maximum u-velocity should be at or near the lid
    println("   Max u-velocity: $(@sprintf("%.3f", u_stats.max)) (should be ~1.0 for lid-driven flow)")
    
    # Velocity magnitude distribution
    speed = sqrt.(result.u.^2 .+ result.v.^2 .+ result.w.^2)
    speed_max = maximum(speed)
    println("   Max speed: $(@sprintf("%.3f", speed_max))")
    
    # Check stream function and vorticity ranges
    if size(result.ψ, 4) == 3  # Vector stream function
        for k in 1:3
            ψ_comp = result.ψ[:,:,:,k]
            ψ_finite = ψ_comp[.!isnan.(ψ_comp)]
            if length(ψ_finite) > 0
                println("   ψ[$k] range: [$(@sprintf("%.3e", minimum(ψ_finite))), $(@sprintf("%.3e", maximum(ψ_finite)))]")
            end
        end
    end
    
    if size(result.ω, 4) == 3  # Vector vorticity
        for k in 1:3
            ω_comp = result.ω[:,:,:,k]
            ω_finite = ω_comp[.!isnan.(ω_comp)]
            if length(ω_finite) > 0
                println("   ω[$k] range: [$(@sprintf("%.3e", minimum(ω_finite))), $(@sprintf("%.3e", maximum(ω_finite)))]")
            end
        end
    end
    
    # Convergence behavior
    if length(result.convergence_history) > 1
        println("\n📈 Convergence:")
        println("   Initial residual: $(@sprintf("%.3e", result.convergence_history[1]))")
        println("   Final residual: $(@sprintf("%.3e", result.convergence_history[end]))")
        
        if length(result.convergence_history) >= 10
            mid_idx = div(length(result.convergence_history), 2)
            println("   Mid-point residual: $(@sprintf("%.3e", result.convergence_history[mid_idx]))")
        end
        
        # Check for monotonic convergence (at least in later iterations)
        if length(result.convergence_history) >= 5
            last_5 = result.convergence_history[end-4:end]
            is_decreasing = all(last_5[i+1] <= last_5[i] * 1.1 for i in 1:4)  # Allow 10% tolerance
            println("   Monotonic convergence (last 5): $(is_decreasing ? "✓" : "✗")")
        end
    end
    
    # Check for NaN or Inf values
    has_nan_u = any(isnan.(result.u))
    has_nan_v = any(isnan.(result.v)) 
    has_nan_w = any(isnan.(result.w))
    has_inf = any(isinf.(result.u)) || any(isinf.(result.v)) || any(isinf.(result.w))
    
    println("\n🔍 Solution health:")
    println("   Contains NaN: $(has_nan_u || has_nan_v || has_nan_w ? "⚠️" : "✅") ($(has_nan_u ? "u " : "")$(has_nan_v ? "v " : "")$(has_nan_w ? "w" : ""))")
    println("   Contains Inf: $(has_inf ? "❌" : "✅")")
    
    # Performance metrics
    if result.solve_time > 0
        dof_estimate = (length(result.x)-2)^3 * 0.6  # Estimate DOF accounting for L-shape
        println("\n⚡ Performance:")
        println("   Time per iteration: $(@sprintf("%.3f", result.solve_time / result.iterations))s")
        println("   Time per DOF: $(@sprintf("%.1e", result.solve_time / dof_estimate))s")
        
        if result.solve_time > 60
            println("   ⚠️  Long solve time - consider smaller problem or optimizations")
        end
    end
end

function main()
    println("3D Lid-Driven Cavity Solver Test")
    println("Testing mathematical formulation and numerical implementation")
    println("Julia version: $(VERSION)")
    
    # Check if we can load the package properly
    try
        lid_driven_info()
    catch e
        println("Warning: Could not display package info: $e")
    end
    
    # Run basic test
    success, result = test_3d_solver_basic()
    
    if success
        println("\n🎉 SUCCESS: 3D solver test completed!")
        
        # Optional: Run additional validation tests
        println("\n💡 Suggestion: For further validation, compare with:")
        println("   - 2D solver results (should see similar flow patterns)")
        println("   - Published lid-driven cavity benchmarks")
        println("   - Python reference implementation")
        
        return 0
    else
        println("\n💥 FAILURE: 3D solver test failed!")
        println("   Check the implementation for errors.")
        println("   Consider debugging with smaller grid sizes.")
        
        return 1
    end
end

# Run the test if this script is executed directly
if abspath(PROGRAM_FILE) == @__FILE__
    exit(main())
end